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Abstract 

The coupled system, where one is a degenerate parabolic equation 
and the other has not a diffusion term arises in the modeling of European 
options with liquidity shocks. Two implicit-explicit (IMEX) schemes that 
preserve the positivity of the differential problem solution are constructed 
and analyzed. Numerical experiments confirm the theoretical results and 
illustrate the high accuracy and efficiency of the schemes in combination 
with Richardson extrapolation 
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1 Introduction 

We study numerically a parabolic-ODE system modelling option pricing liquid¬ 
ity shocks . The presence of liquidity shocks is a source of non-liquidity risk and 
makes this market incomplete. Ludkowsky and Shen [5] investigate a nonlinear 
pricing mechanism based on utility maximization. They consider the investor 
whose utility is described by an exponential utility function 

U{x) = -e-^b (1) 

where 7 > 0 is the coefficient of risk aversion. The investor seeks to maximise 
utility of both terminal wealth and option payoff at time horisont T < 00 , which 
is chosen to coincide with the expiration date of all securities in market model. 
Properties of the exponential utility function (1) imply that the value functions 
can be expressed as 

U\t,X,S) = = 0,1, (2) 

where A = Xt is the wealth process and the functions i?*(b S) are related to the 
price of options in the two states, see (7) below. Then the pair {R^{t,S), i = 
0 , 1 } is the unique viscosity solutions of the coupled semi-linear system, 


_ £10 I _ Q 

t y y ' 


= 0 , 


( 3 ) 


1 



The terminal conditions are: 


R\T,S) = h{S), 1 = 0,1. (4) 

Here cr is volatility of the underlying, are transition intensities from 

state ( 0 ) to state ( 1 ) and vice versa, respectively, fj, is drift of the underlying 
and do = see [5] for more details. 

Using [/* and U* , the buyer’s indifference price p (initial state 0) and q 
(initial state 1 ) are defined via 

U°it,X-p,S) = V%t,X), U\t,X-q,S) = V\t,X), (5) 

where U, V are the optimal solutions for terminal wealth with and without 
options respectively. The value functions U*, i = 0,1 are given by 

U* = e-T'^Ui(t) and U*(t,X,5') = i = 1,2 ( 6 ) 

and the functions Fo{t),Fi{t) by, 


Fo{t) = + C 2 e^'^\ 

Fi{t) = —{ci((io + ^^01 ~ + C2{do + t'oi — A 2 )e^’^*, 

r'oi 


where 


Ai,2 — 


do + r'oi + ^^10 i a /{do + r'oi + ~ 


A 2 do -\,T 1 Ai do -\^T 

Cl = --—e 1 and C 2 = -r- ^ . 

A 2 ~ Ai Ai — A 2 

Then, we obtain from (2), (5), ( 6 ) 

p = i?o+7-ilnUo(t), g = i?^+7-MnFi(<) 
and from (3), (4) the parabolic-ordinary system for p and q 


Pt 

qt 


+ h^'^S^pss - = 0 , 


_ 'Hio. Fo_^—'l{q—p) _j_ 210. 1 ^1 _ Q 


7 Fo 


7 Fi 


7 ^^1 


with terminal conditions 


(7) 

( 8 ) 


p{T,S)=q{T,S) = h{S). (9) 

The numerical solution of the system ( 8 ) is the main object of the present 
paper. The numerical treatment of the boundary layer effect for small values of 
CTo and 7 , the degeneracy at S' = 0 of the parabolic equation and the exponential 
nonlinearity lead to challenging problems [10]. The introduction of exponential 
nonlinear terms is an available assumption based on the financial nature of 
the model system ( 8 ). There are many numerical schemes to solve nonlinear 
parabolic and hyperbolic equations. However, very few have dealt with an 
exponential nonlinear term. The special nature of the nonlinear exponential 
term for a hyperbolic problem is discussed in [10]. A possible way to build an 
efhcient numerical solution of ( 8 ), (9) is to implement an IMEX method [1,9] 
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. In this procedure the diffusion term is discretized implicitly in time and the 
reaction terms are discretized explicitly. 

An IMEX method for numerical solution of reaction-diffusion equation with 
pure Neumann boundary conditions is developed in [3]. IMEX schemes, by 
applying explicit approximation both integral term and the convection term 
and an implicit approximation for the second differential term are developed for 
integro-differential equations of finance in [2]. 

The rest of the paper is organized as follows. In the next section some results 
concerning well-posedness of Cauchy problem for Q and a comparison principle, 
obtained in [4] , are discussed . Also, two lemmas concerning discrete maximum 
principle [6,7] are formulated. In Section 3 an implicit-explicit linear scheme is 
introduced. Comparison discrete principle and convergence of the scheme are 
proved. Similar results are obtained for an IMEX linearized scheme in Section 
4. The computational experiments in Section 5 confirm the applicability of 
our schemes and the theoretical results. Finally, Section 6 summarizes our 
conclusions. 

Notation Let be a bounded interval in = (0, oo) and let 
denote the space of continuous functions on H with the norm of any w G 
defined by ||ry||n = sup,j,gQ |ic(a;)|. For each integer fc > 1 , let (7^(11) denote the 
space of fc - times differentiable functions on with continuous derivatives up 
to and including those of order k, with the norm of any w G (7^(12) defined by 
||w||fc,n = maxo<i<fc ||o. The notational conventions |w|o,n = ||w||o,n = 

||w^||n = 11^11 are adopted . The explicit reference to fl is dropped whenever the 
domain question is evident. For any mesh functions on arbitrary mesh = 
n = the discrete maximum norm is defined by = 

maxo<j<Ar |wi|. 

Maximum norms and semi-norms for smooth functions of two variables 
are introduced in a similar way. Let Qt = (0,T) x fl. Then ||ru||Q.j, = 
supj-j, |u>(a;,r)| and if C^{Qt) is the space of all functions on Qt with 
continuous derivatives then 

c (Qt) = g C°{Qt) for i, j = 0,1, 2 ,... with 0<i + 2j<k 

2 Preliminaries 

In this section we will describe some properties of the solution to system (8) 
using results obtained in [4]. Also , following [6,7], two lemmas , concerning 
discrete maximum principle (DM) are formulated. 

We will consider solutions of (8) satisfying 

\p \) kl) \h\ < A exp (aln^ S) = (10) 

for some positive constants A and a. In [4], well-posedness in weighted Sobolev 
spaces and comparison principle for the corresponding Cauchy problem (8), (9) 
are established. With sufficient smoothness of the initial data the weak solutions 
are classical ones. 

In this paper we use the comparison principle for classical solutions p{S, t), q{S, t) 
of the problem (8), (9), i.e p G <7((0, -l-oo) x (o, r])n(7^7((0, -foo) x (0, T)), q G 
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C'((0,+(^)^x (0,T]), qt € ((0,+oo) X (0,r)). 

Proposition 1 ([4]) Let {pi,qi) and {po, qo) be two classical solution of problem 
(8), (9) corresponding to terminal data h = hi{S) and h = ho{S) , respectively. 
If there exists some positive constants A and a such that pi{S,t) and hi{S), i = 
0,1 satisfy conditions (10), then 

- ho) < pi{S,t) -po{S,t) < sup(ft,i - ho), 
mi{hi-ho) <qiiS,t)-qo{S,t) <sup{hi-ho), 

In particular, let h{S) be hounded from below (or from above) by a constant 
h{S) > hi, (resp. h{S) < h* and the pair p(S,t),q{S,t) he a classical solution 
of the terminal problem (8), (9). Then 

piSjt) > hi, and q{S,t) > hi, (respectively p{S,t) < h* and q{S,t) < h*). 
for any S € (0, +oo) and any t S (0, T]. 

By making the substitutions t = T — t, u = yi?® 

(3) becomes 

LP{u, v) =Ut — ^a'^S'^uss + ae“e“'^ 

LP{u,v) =Vr + ce’'e““ — c = 0, 

where a = voi, b = do + voi, c = J^io- In accordance with (9) we take the initial 
conditions to be 

u{0,S) = uo{S)=jh{S), v{0,S)=vo{S)=jh{S). (13) 

For a call option, 


and V = the system 


-5 = 0, 


( 12 ) 


h{S) = ina,x{S — K, 0). (14) 

We assume ground conditions for u, v of the form (10). In the next sections, 
the analysis of the difference approximations of problem(12)-(14) will use the 
following comparison principle that follows from those one for (p, q): 
Proposition 2 Let {u,v), {u,v) G C([0, T) x (0, +oo)) fl (^^’^((O, T) x (0, +oo)) 
be two pairs of classical solutions of ( 12 )-( 14 ) corresponding to the initial data 
h = h and h = h, respectively and such that conditions of the type (10) hold. If 
the following inequalities also hold: 

L^{u,v) > L^{u,v), L°(u,v) > L^{u,v) and h > h, (15) 

then 

u > u,v > V. 

Hereinbelow we will use the following canonical form of writing a 3-point 
difference scheme 


Aiyi—i CiPi “t“ — Fi, i — 1 , 2 ,..., N 1 

2/0 = hi, VN = 02- 


(16) 


The discrete comparison principle for problem (16) was proved in [6,7] and 
is formulated in the following way. 
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Lemma 2.1. Let the conditions 


A>0, B,>0, D, = Ci-A,-B,>0, i = l,2,...,N-l (17) 

be fulfilled. Then the solution of the difference scheme (15) satisfies the inequal¬ 
ities 


> 0, i = 0,...,N, if Fi>Q, i = l,...,N -I, /xi > 0, ^^2 > 0; 


2 /i < 0, i = 0,...,N, ifFi<0, i = l,...,N-l, /xi < 0, ^2 < 0. 
Lemma 2.2. Let the conditions 


|A,|>0, |B,|>0, A = |Ci| - 1^,1 - |B,| > 0, i = l,...,N-l 
be met. Then for the solution the problem (16) the estimate holds 


|| 2 /|lc(n'^) < max<^ IM1I.IM2I 


c(a'V) 


3 Implicit-Explicit Linear Scheme 

In this section, we develop a linear IMEX scheme to solve the coupled semi-linear 
parabolic-ordinary system problem (11)-(12). 

For call option one possible pair of boundary conditions is, see e.g. [10,11] 

u(r, 0) = (pfir) = 0, u{t, S) = Prij) « S'max for large S. (18) 
The left natural boundary condition for u is 

Ur{T, 0) = -b b. (19) 

On Qt = O X [0, T] we introduce the uniform mesh wst = ws x Wt '■ 


ws = {Si=iAS, AS > 0, 1 = 0,1 ,...,/; IAS = Smax}, wg = wsU{5'o,S'/}; 

Wr = {Tj=jAT, Ar > 0, j = 0,1,...,J; JAr = r}., Wr = Wr U {to,tj} 

On the discrete domain wst we approximate the problem (12)-(14) by the 
difference scheme 


LP{U,V) = 


A+i _ Tji 1 — 2A^^ -b 

Oi Oj i 2r.2^i-l ^ i+1 


At 


^ A c 2 

- 2 ^ S 


{Asy 


ae 


-vlA 


-6 = 0 , 
( 20 ) 
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y'i+1 _ 

L^(U,V) = ^ - c = 0, i = 0,l,...,/, (21) 

l\T 

j = 0,1,..J- 1; 


U° = Uo{Si), i = 0,l,...,/; 


( 22 ) 


U^=(pi{Tj), Uj=ipr{Tj), j = 0,(23) 

V° = Vo{Si), i = (24) 

The natural boundary condition can be approximated as follows 

- AT(ae“e-^o - 5). 

On the (j + l)-th, j = 0,1) • ■ •, — 1 time level the scheme (20)-(23) has the 

form 


-AdJt^ + edit" - BdJitl = Fi, 

= vi - Arce-^'+^i + c, 


where 


1 


Ai = Bi = 


2 (AS') 


e2 1 

1 Ci = -7 - Ai + Bi, 


At 


F, = —U^-ae-^'e^'+b, i = 1,..., I - 1; 

lAT 

For the truncation error corresponding to (20) we find 

1 3‘^u 

Trl = -At—(t,-+i - 6»i At, S,) 


(25) 


+At ( ^(a +1 “ P AT,Si_i) + ^(a +1 “ pAr^Si) + ^(a+i “ p’''At, S^+i) 


- 0i-AS)| 

+At—(t,+i -pAT,S,)eA-^+i.«^) = _(r^.+i _ pAt, Si)e-A-^+i.sd 

= 0(At) + (AS)2. 
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For the truncation error corresponding to (21) we get 


1 3^1! 3ii 

Tr2 = -Ar^(rj-+i - 6 » 2 Ar, S,) + At—(tj+i - rjAr, 

- —(r,+i - = 0 (At). 

0 < p, p“,p+,p, P < 1, 0 < 01 , 9 ^, 9 ^ <1. 

In accordance with Notation we define the strong norms on the meshes Ws 
and wst, respectively, 

ll^llc(Trs) = max |2i|, ||z||c(T7s.) = max \zf\. 

0 <i<l o<t<i 

- 0 <i<J 

Let denote 


Cu= sup |m(t,S')|, Cy = sup |u(t,S')|. 

(r,S)GQT (T,S)eQT 

Theorem 1 Suppose that there exists classical solution {u,v) € C^'^(Qt) of 
problem (10)-(If)- Then for sufficiently small AS and At the following error 
estimate holds: 


Ik “ U\\c{ws.r) + Ik “ V\\c{ws^) ^ C{At + {ASY), 

where the constant C doesn’t depend of AS and At. 

Proof Define the errors sj, pf by 

el=ul -u{Tj,Si), pl=Vf -v{Tj,Si), i = l,...,I. 
Then {s(}, {p)} satisfy the linear system of algebraic equations: 

=Fi, i = 1,1, 


(26) 


A+i — 0 A+i — Q 


where 


and 


1 

At 


F) = —— s) + a] 


+ At PI- 

Here a) and P) are the local truncation errors corresponding to the difference 
equations (20) and (21), respectively. They will be estimated as follows. 

Let us derive the truncation error corresponding to nonlinear (right) part: 
For the nonlinear, right hand side of the first equation we obtain 
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= {l-^4+ Oiinim +si+ 


= (l+ei- - 0(£^/i^) + Odeif) + Od^f), 

and another one form 

e-v/e'^i = + 0(ej) + 0(^j). 

Now, taking into account the Trl, we have 

In a similar way, we find 


/3f = 0 (At) + (fd - ef)e-“("^- 
Applying Lemma 2.1 we get 


+ 0(sifd)+Odsif) + OdMif). 


llei+dic < AtIIFII, 

where || • ||c is the strong norms C{ws) as defined above. 

We estimate HA-’+^H : 

+0(Ar) + OdASf) + OiWs^r) + 0(1^^)). 

Next, 


\\^^+d\<y^\+^rie^-^d\ed\+OiAr)+0{\\s^r)+0{\\^ 

Therefore, 

||A+i|| + ||M^'+'||<(l + 2ATe^“e^-)||5^|| 


+ (1 + ArCMWfi^W + AT(0(Ar) + (AS)^ + 0{\\A\d) + OdlM'lH)- 
For j = 0 we have £° = 0, Mi = 0 and then 

a° = 0(Ar) + (AA)2, /jO = 0 (At). 

Since ||£°|| = ||m°|| = 0 , we get 

||£d| = CAt{At + {ASf), ||£d| = C'Ar(AT + (AS')^). 
Therefore, by induction we have 



+ iim^'+'ii < (1 + 2Arc„a)(ik^ii + iim^ii) 


+ArC'(Ar + (AS')2). 

which implies that 

||A+1|| + y^+^\\ < C^(l + 2AtCM^At(At + (ASf) 

k^O 

C(At + (A5)2)(((1 + 2ATC7,a)'^(l + 2ArC,a)-i 
< C'(At + (AS')2)n 

The following discrete comparison principle for the {U, V) is crucial for the 
positivity of the discrete approximations of the indeference prices p and q on 
the base of the scheme (20)-(24). 

Theorem 2 Let the assumptions of Theorem 1 hold Let also {U,V) , {U_,Y.) be 


grid functions defined on wst and the inequalities hold: 

LP(U,V) > LP{U,V), L°(U,V) > L°{U,V), (17) 

V°>Vl 1 = 0,...,1, (28) 

vi>Vl Um>U%, J = (29) 

Then for sufficiently small AS and At we have 

Ul>m^ Vl>Yf, i = 0,l...,I, j = 0,l,...,J. (30) 


Proof. Let introduce 


yl = ul-m. 4=yl-Yi, i = 0,i...,i, J = 0,1,...,J. 
Then, from (26) we obtain 


yj -yj 

At 


- 




-^yi(0)>O, (31) 




-|-c(e ^*e^*—e —^e—^)>0, i = l..j = l,...,J—1, (32) 


Using the mean-value theorem we get 

e-yle^i -e-yie^i =e-^'e^'(jjl-zi), 

ui =m + Wl-m), ^ = vl+e(vl-vi),o<e<i, 

ui=m+mi-m),y^ = yl+o(yi-vi), o<e<i. 
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We rewrite (30) in the form 


Aytl - ayr^B.ylX^ > -F„ 


,J + 1 R.^,t + l 


Fi = -^yi - “e (yi -^D- 


Next, we rewrite (31) in the form 
.J'+i / 1 


At 


> 


At 


-u^ v^ 
— ce '6' 


z{+ce yi. 


(33) 


We apply the method of mathematical induction with respect to j to prove 


that 


yf >0, 4>0, t = 0,l,...,/, j = 0,l,...,J. (34) 

From (16), (17), we have 

y° >0, 2° > 0, i = 0,1 ,..., /, 

Assuming that (32) holds when j = k — \, we will show that for j = k the above 
inequalities are true. 

On the base of Theorem 1 we can confirm that for sufficiently small At, AS 
there exists constants Cu,Cv, such that 

max(||[/||, \\U) < 2Cu, max(||y||, ||£) < 2C'„. 

Then, if it necessary, we choose At in additional smaller such that 

At < min(a, c)e^^“e^‘"” (35) 

By induction, y^~^ > 0, > 0 and using (Lemma 2.1) we conclude that 

Fi > 0, i = 0,1,...,/— 1. Now Lemma 2.1 implies yf > 0, i = 0,1,..., /. It is 
clear from (32) and (34) that zf > 0, i = 1,..., 7 — 1. □ 


4 Implicit-Explicit Linearised Scheme 

Let us consider first the implicit scheme: 

--2^ -h ae • e • -6 = 0, 


{ASf 


yr - yi , 

—— - -+ce - c = 0, 

At 

z= 1,2,...,/-!; j = 0,l,..., J-1. 


with boundary and initial approximations (22)-(24). 

By Taylor expansion we get 

= e-vf gc/f (1 + ^ + g-WgC/f ^,+1 _ ^+1) 


(36) 

(37) 
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+o((c/|+i - ulf ) + o((y/+i - V^f), 

= e-^^e^(l -V^ + Ul) + +1 - C/|+i) 

+o((c//+i - Ulf) + o((y/+i - V^f). 

We drop the O-terms and the results we insert in (36) and (37) to obtain: 


-AdJtl + C^Ul+^ - mlfl + = b, (38) 

= Gi, (39) 


where 


Ai — Bi — -a^ ~ Ar ’ 

4 = -are^i-vf p, = -^Ul - are^^-^'(l + 1^^' - U/) + &At, 
At 

= -cAre^^'Ki = 

At 

G. = -^Vf - cre^-^^ (i _ + C/|) + c. 


Since > 0 , the diagonal domination can significally increase in 

comparison with IMEX linear scheme, see system (20),(21). 

Theorem 3 Let the assumptions of Theorem 1 hold. Then suppose that there 
exists classical solution {u,v) € C^ '^{Qt) of problem (10). Then for sufficiently 
small AS and At the following error estimate holds: 


l|u - UWciwsP + Ik - ^llc(t«sA < C;{At + {AS)b, 

where the constant C doesn’t depend of AS and At. 

Proof. Substituting from (37) into (36) the first one we get 


-AiUlfl+ Ci- 


DiEi 

K,^ 


— F- — -F- 


If- 


Ki ki " ' ’ ’ 


with Uf f = 0,1,..., / and Uf Uf j = 0,1,..., J given by (22),(23) and (24). 
For the error we have the linear system of algebraic equations 


-AiSlfl + ( G, - 


DiEi 


K, 


i_A \ B+i 


~ k+1 — + CC 


A+1 


= 0, 4+^ = 0 
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i+1 

R - 


= e: 


i+i 


K,, 


A 

At 


+ Pa * = 1 


,/-l. 


Further we follow the line of Theorem 1 to complete the proof. □ 

The scheme (36), (37) also has similar comparison properties of the linear 
IMEX scheme described in Theorem 2. 


5 Numerical Experiments 

In the section we perform numerical experiments to illustrate the accuracy, effec¬ 
tiveness and convergence of the implicit-explicit linear scheme (20)-(24) (Scheme 
1) and implicit-explicit linearized scheme (38),(39) (Scheme 2) developed in this 
article. We provide experiments both with uniform and non-uniform meshes. 
Also, we present results of numerical experiments using Richardson extrapola¬ 
tion in time. 

The Tables (presented results) show the accuracy in maximal discrete norm 
II • II and convergence rate at final time T, using two consecutive meshes with 
formulas 


Ratio = log2(E7/2/E/), = \\We:o - Wiloo, 

where Wex and W are the exact and the corresponding numerical solutions, 
respectively. In our case Wex is or . 

In Tables [TJ [2l we give the results from the computations IMEX linear 
Scheme I. 


Table 1: Convergence results for at the money {S = 2,K = 2,5'„iin = 0 and 
Smax = 5) and At = AS/2 based on Scheme 1 


/ 

Value 

Difference 

Ratio 

Value 

R^ 

Difference 

Ratio 

30 

0.246669 



0.235165 



60 

0.247438 

7.70e-04 


0.235917 

7.52e-04 


120 

0.247749 

3.11e-04 

2.48 (1.31) 

0.236218 

3.01e-04 

2.50 (1.32) 

240 

0.247887 

1.38e-04 

2.25 (1.17) 

0.236349 

1.31e-04 

2.30 (1.20) 

480 

0.247952 

6.50e-05 

2.12 (1.09) 

0.236410 

6.10e-05 

2.15 (1.10) 

960 

0.247983 

3.10e-05 

2.10 (1.07) 

0.236439 

2.90e-05 

2.10 (1.07) 


Table 2 is based on a non-uniform grid and also shows that the scheme is 
first order in time. Here we use Tavella-Randal [8] mesh: 

Si = K + a (^C 2 j + Cl , 

. ,-l fSxnin-K\ . ,.i fSxaax-K\ 

Cl = smh - I C2 = smh - . 

\ a J \ a J 

In this case, we choose to concentrate mesh points around the strike price K 
since we expect the error to be largest there. In Table 3 we list the results from 
computation with Scheme 2 that for this non-uniform grid the results are still 
first order accurate in time as in the uniform case. 
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Table 2: Convergence results for at the money {S = 2, K — 2, Snun = 0 and 
i^max = 5) and taking Ar = ASi/2 and using nonuniform Tavella-Randal grid 
with Q = 15 based on Scheme 1 


I 

Value 

Difference 

Ratio 

Value 

Difference 

Ratio 

30 

0.247196 



0.235660 



60 

0.247863 

6.67e-04 


0.236305 

6.45e-04 


120 

0.248124 

2.61e-04 

2.56 (1.35) 

0.236552 

2.47e-04 

2.61 (1.38) 

240 

0.248238 

1.14e-04 

2.29 (1.20) 

0.236657 

1.05e-04 

2.35 (1.23) 

480 

0.248291 

5.30e-05 

2.15 (1.10) 

0.236706 

4.90e-05 

2.14 (1.10) 

960 

0.248322 

3.10e-05 

1.71 (0.77) 

0.236735 

2.90e-05 

1.69 (0.76) 


Table 3: Convergence results for at the money {S = 2, K = 2, ^min = 0 and 
Smax = 5) and taking Ar = AS/2 based on Scheme 2 


I 

Value 

R^ 

Difference 

Ratio 

Value 

R^ 

Difference 

Ratio 

30 

0.246685 



0.234952 



60 

0.247444 

7.59e-04 


0.235812 

8.60e-04 


120 

0.247752 

3.08e-04 

2.46 (1.30) 

0.236165 

3.53e-04 

2.44 (1.28) 

240 

0.247889 

1.37e-04 

2.25 (1.17) 

0.236323 

1.58e-04 

2.23 (1.16) 

480 

0.247953 

6.40e-05 

2.14 (1.10) 

0.236397 

7.40e-05 

2.14 (1.09) 

960 

0.247984 

3.10e-05 

2.06 (1.05) 

0.236433 

3.60e-05 

2.06 (1.04) 


Now, we improve the convergence in time applying Richardson extrapolation 
[4]. To this aim we use the formula 

_ 2PWn - Zn 
” 2P - 1 

where p is order of numerical solution (1 in our case) and Wn is the solution 
obtained using time step At/2 and is the solution obtained using time step 
At. The resulting solution has order of accuracy p + 1 [4]. Table [5] shows 
the result of applying this technique to the Scheme 1. The order of accuracy 
in time is now two. Similarly this technique is applied to Scheme 2, see Table 
6 . Hence the convergence is much slower but smoother compared to the explicit 
based Seheme 1 due the error of linearisation. The tables shows second order 
in time. 
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Table 4: Convergence results for at the money {S = 2, K — 2, Snun = 0 and 
i^max = 5) and taking Ar = ASi/2 and using nonuniform Tavella-Randal grid 
with Q = 15 based on scheme 2 


I 

Value 

Difference 

Ratio 

Value 

Difference 

Ratio 

30 

0.248722 



0.237005 



60 

0.249432 

7.10e-04 


0.237812 

8.07e-04 


120 

0.249715 

2.83e-04 

2.51 (1.33) 

0.238139 

3.27e-04 

2.47 (1.30) 

240 

0.249839 

1.24e-04 

2.28 (1.19) 

0.238283 

1.44e-04 

2.27 (1.18) 

480 

0.249897 

5.80e-05 

2.14 (1.10) 

0.238351 

6.80e-05 

2.12 (1.08) 

960 

0.249928 

3.10e-05 

1.87 (0.90) 

0.238387 

3.60e-05 

1.89 (0.92) 


Table 5: Convergence results for at the money {S = 2, K = 2, ^min = 0 and 
‘S'max = 5) and taking At = AS/2 based on Scheme 1 using Richardson extrap¬ 
olation _ 


I 

Zn, Wn 

v„ 

Difference 

Ratio (order) 

10 

0.2451080 




20 

0.2465578 

0.2480075 



40 

0.2472811 

0.2480045 

3.02e-6 


80 

0.2476431 

0.2480051 

5.79e-7 

5.22 (2.38) 

160 

0.2478242 

0.2480053 

1.96e-7 

2.96 (1.56) 

320 

0.2479148 

0.2480053 

5.13e-8 

3.82 (1.93) 

640 

0.2479600 

0.2480053 

1.27e-8 

4.05 (2.02) 

1280 

0.2479827 

0.2480053 

3.08e-9 

4.12 (2.04) 

2560 

0.2479940 

0.2480053 

7.45e-10 

4.13 (2.05) 


In Figure [T] we compare options values p and q at issue and maturity in the 
liquid and illiquid states using the parameters p = 0.06, tT=0.3, r'oi = 1, vio = 
12, RT = 2, T = 1, 5'max = 5 and 7 = 1 using the Scheme 1 . Figure [5] illustrate 
the linearised scheme, using the same parameters. Figures 1,2 illustrate the 
positivity of the solution {p,q), using both schemes. 


6 Conclusions 

In this work we have considered one-dimensional problem of European options with liq¬ 
uidity shocks. We have constructed and analyzed two IMEX finite difference schemes 
that preserve the positivity property of the differential solution. The second one(the 
IMEX linearized scheme) has better diagonal domination, respectively monotonicity. 
It would be interesting to consider extensions of the IMEX schemes to the American 
options with liquidity shocks. In this case one has to solve a free boundary problem. It 
could be written as a linear complementary problem which could be discretized using 
the schemes given here. The extension is beyond the scope of this paper, and we leave 
it for further work. 
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Table 6: Convergence results for at the money {S = 2, K = 2, S^in = 0 and 
Smax = 5) and taking At = AS/2 based on Scheme 2 using Richardson extrap¬ 
olation. _ 


I 

Zn, Wn 

Yn 

Difference 

Ratio (order) 

10 

0.2451717 




20 

0.2465832 

0.2479947 



40 

0.2472928 

0.2480023 

7.64e-6 


80 

0.2476486 

0.2480045 

2.14e-6 

3.57 (1.84) 

160 

0.2478269 

0.2480051 

6.22e-7 

3.44 (1.78) 

320 

0.2479161 

0.2480053 

1.78e-7 

3.49 (1.81) 

640 

0.2479607 

0.2480053 

4.93e-8 

3.62 (1.85) 

1280 

0.2479830 

0.2480053 

1.32e-8 

3.73 (1.90) 

2560 

0.2479942 

0.2480053 

3.46e-9 

3.81 (1.93) 

5120 

0.2479998 

0.2480053 

8.95e-10 

3.87 (1.95) 

10240 

0.2480026 

0.2480053 

2.29e-10 

3.91 (1.97) 




(a) p at t = 0 and t = T (b) q at t = 0 and t = T 


Figure 1: Comparing European option values at issue and maturity in the liquid 
and illiquid states for the IMEX Linear scheme 
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Figure 2: Comparing European option values at issue and maturity in the liquid 
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